function my_map2draw(lat,lon,raster,z,option)
% my_map2draw(lat,lon,image,z,option)
% -------- 自定义函数 --------
% 带经纬度绘制raster，使用pcolor
% lat：纬度向量
% lon：经度向量
% image: 是2维空间矩阵，现支持第三维是时间的三维矩阵，上北下南左西右东，
%   需要行数和lat大小一致，列数和lon大小一致。
% z：    可选参数，double，第三维的下标。
% option:通过键值对模式输入，如func("KeyName",KeyValue)
%   "Center": true/false, Climit要中心化，例如[-1,1]这样对称，能够凸显出相对0的差距
%   "CenterType": 1/2, 1是按照绝对值最小值中心化，默认，2是按照绝对值最大值中心化
%   "Climit": [cmin,cmax] 设定caxis()的climit
%   "SetNaN": true/false 是否需要调用SetNaN函数
%   "NaNTag": 需要设为NaN的值，默认为0
%   "ColorMap":string, colorbrewer提供的颜色主题名称，通过在主题名称前面加*
%       使得颜色图翻转。默认颜色图是"*Spectral"。
%   "ColorDisperse": true/false,颜色是否离散，默认连续
%
% 函数方法说明：
%   使用pcolor达到和image相同的绘制效果。选用pcolor的原因是他不绘制nan数据，
%       而imagesc会讲nan数据也给定颜色。
%   pcolor绘制的时候是以左下为坐标起点的，使得绘制结果和image上下颠倒，
%       用flipud进行矩阵的倒置。
%   

% 作者：sid(cnqdcyq@163.com)
% log：
%   2021/11/20 ver01 create，根据map2draw改编而来
%              ver1.1 update: 增加了SetNaN的功能
%              ver1.2 update: 增加了raster必须是数值类型的要求,增加对全NaN值
%                   图的支持
%   2021/11/21 ver1.3 update: 增加了对颜色图brewermap函数的支持。完善了
%                   help信息。
%	2022/01/05 ver1.4 fixbug: 增加了再连续变量时，Center的支持功能
%	2022/01/07 ver1.5 update: 增加了CenterType的功能
%	2022/01/08 ver1.6 update: 修改了经纬度的格式，带度，去掉经纬度名称

	arguments
		lat
		lon
		raster {mustBeNumeric,mustMatchCoor(lat,lon,raster)}
		z = 1
		option.Climit = []
		option.Center = false
		option.CenterType double = 1
		option.SetNaN = false
		option.NaNTag = 0
		option.ColorMap string = "*Spectral"
		option.ColorDisperse {mustBeNumericOrLogical} = false 
	end

    % if SetNan is true, set NaNTag in matrix to NaN
    if option.SetNaN
        raster = my_setNaN(raster,option.NaNTag);
    end
    
    if option.ColorDisperse
        SetValue = my_getunique(raster(:,:,z),'omitnan');
        Climit = getClimit(option.Climit,SetValue);
        Type = returnType(Climit,SetValue);
        Climit = getNewClimit(Climit,Type,...
			option.Center,option.CenterType);
        cmap = getCmap(Climit,option.ColorMap);
    else
        Climit = option.Climit;
        cmap = brewermap([],char(option.ColorMap));
    end
    colormap(cmap)

    % plot spatial pattern
    pcolor(lon,flip(lat),flipud(raster(:,:,z)));
    
	% change format of axis to degrees
	xtickformat('degrees');
	ytickformat('degrees');
    shading flat
    axis equal
    
    if ~isempty(Climit)
		% climit not empty means colordisperse,set caxis
        caxis(Climit);
	end


	if ~option.ColorDisperse 
        % if color is continous
		if ~isempty(Climit) % get Climit
			% if Climit has been set
			caxis(Climit);
		else
			% if Climit not set
			Climit = caxis;
		end

		if option.Center && Climit(1) * Climit(2) < 1 
            % if Center 
			Climit = getCentered(Climit,option.CenterType);
        end
		% set caxis
		caxis(Climit);
    end
    
    if ~option.ColorDisperse || numel(SetValue) > 1
        colorbar;
    end
end

function mustMatchCoor(lat,lon,raster)
if any(size(raster,1,2)~=[numel(lat),numel(lon)])
    error('经纬度向量的大小和矩阵的大小不匹配！')
end
end

function output = getClimit(Climit,SetValue)
    oldCaxis = [min(SetValue),max(SetValue)];
    if isempty(Climit)
        output = oldCaxis;
    else
        output = Climit;
    end     
end

function output = returnType(Climit,SetValue)
if isempty(SetValue)
    output = [0,0];
else
    Tag = SetValue < Climit(1) | SetValue > Climit(2);
    output = [Tag(1),Tag(end)];
end
end

function output = getNewClimit(Climit,Type,Center,CenterType)
    output(1) = Climit(1) - double(Type(1)) - 0.5;
    output(2) = Climit(2) + double(Type(2)) + 0.5;

    if Center && (output(1)*output(2)<0)
		output = getCentered(output,CenterType);
    end
end

function output = getCmap(Climit,ColorMap)
    ClassNumber = Climit(2)-Climit(1);
    output = brewermap(ClassNumber,char(ColorMap));
end

function out = getCentered(Climit,CenterType)
	% get Centered Climit according to CenterType
	switch CenterType
		case 1
			Abs = min(abs(Climit));
			out = [-Abs,+Abs];
		case 2
			Abs = max(abs(Climit));
			out = [-Abs,+Abs];
	end
end
